# Code to replicate
#   Job Insecurity, Economic Resources, and Democratic Backsliding:
#   Evidence from South Korea
#
# Original analyses conducted in R version 4.3.1
#
# Below are the URLs for the data required to replicate the study
#
# https://v-dem.net/data/the-v-dem-dataset/
# https://www.worldvaluessurvey.org/WVSDocumentationWV7.jsp
# https://www.worldvaluessurvey.org/WVSDocumentationWVL.jsp
# https://kgss.skku.edu/kgss/data.do


#####
# V-Dem
# Figure 1

f1 <- readRDS("~/Desktop/Replication/V-Dem-CY-Full+Others-v13.rds")

library(tidyverse)
select <- dplyr::select
f1 <- f1 %>% filter(country_text_id == "KOR", year %in% c(1948:2022)) %>%
  select(year, v2x_polyarchy, v2x_libdem, v2cacamps_mean) %>%
  mutate(v2cacamps_mean = v2cacamps_mean/4) %>% 
  rename("Electoral Democracy" = v2x_polyarchy,
         "Liberal Democracy" = v2x_libdem,
         "Political Polarization" = v2cacamps_mean) %>% 
  pivot_longer(cols = -year, names_to = "index", values_to = "value") %>%
  mutate(index = factor(index, levels = c("Electoral Democracy",
                                          "Liberal Democracy",
                                          "Political Polarization")))

library(ggplot2)
Fig1 <- ggplot(f1) +
  geom_line(aes(x = year, y = value, color = index), lwd = .65) +
  geom_point(aes(x = year, y = value, shape = index), size = 1.9, color = "white") +
  geom_point(aes(x = year, y = value, shape = index, color = index, size = index)) +
  scale_shape_manual(values = c(20, 18, 15)) +
  scale_color_manual(values = c("#888B8D", "#53565A", "#A7A8AA")) +
  scale_size_manual(values = c(1.5, 1.5, 0.9)) +
  scale_x_continuous(breaks = c(seq(1950, 2020, by = 5)), expand = c(.015, .015)) +
  scale_y_continuous(breaks = c(seq(0, 1, by = 0.2))) +
  labs(x = NULL, y = NULL) +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  theme(axis.text.x = element_text(size = 11, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 11, face = 'bold', color = "black")) +
  theme(legend.position = c(.845, .145), legend.title = element_blank(),
        legend.text = element_text(size = 11, face = 'bold'),
        legend.key = element_rect(fill = "transparent"),
        legend.background = element_rect(fill = 'transparent'))

Fig1


#####
# WVS
# Figure 2, Figure 3, Figure 4
# Table A5, Table A7, Table A8, Table A10, Table A11, Table A12, Table A14, Table A15
# Figure A1, Figure A2, Figure A3, Figure A4, Figure A6

library(haven)
wvs <- read_dta("~/Desktop/Replication/WVS_Wave_7_South_Korea_Stata_v5.0.dta")

library(dplyr)
wvs <- wvs %>% mutate(strldr = 5 - Q235,
                      strl2 = ifelse(strldr < 3, 0, 1),
                      jinsc = 5 - Q142,
                      KR_A24 = ifelse(KR_A24 < 1, NA, KR_A24),
                      outs = ifelse(KR_A24 == 2, 1, ifelse(KR_A24 == 1, 0, NA)),
                      aset = KR_A30,
                      inc = Q288,
                      age = X003R,
                      gen = ifelse(Q260 == 2, 0, Q260),
                      edu = Q275,
                      ideo = Q240,
                      fmj = KR_B12, fjh = KR_B13,
                      apolr = abs(fmj - fjh),
                      region = ifelse(N_REGION_WVS == 410001, "Seoul",
                                      ifelse(N_REGION_WVS == 410002, "Busan",
                                             ifelse(N_REGION_WVS == 410003, "Daegu",
                                                    ifelse(N_REGION_WVS == 410004, "Incheon",
                                                           ifelse(N_REGION_WVS == 410005, "Gwangju",
                                                                  ifelse(N_REGION_WVS == 410006, "Daejeon",
                                                                         ifelse(N_REGION_WVS == 410007, "Gyeonggi",
                                                                                ifelse(N_REGION_WVS == 410008, "Gangwon",
                                                                                       ifelse(N_REGION_WVS %in% c(410009, 410010), "Chungcheong",
                                                                                              ifelse(N_REGION_WVS %in% c(410011, 410012), "Jeolla",
                                                                                                     ifelse(N_REGION_WVS %in% c(410013, 410014), "Gyeongsang",
                                                                                                            ifelse(N_REGION_WVS == 410020, "Ulsan", NA)))))))))))),
                      hap = 5 - Q46,
                      satfin = Q50,
                      tgov = 5 - Q71,
                      tcser = 5 - Q74,
                      tpart = 5 - Q72,
                      tparl = 5 - Q73,
                      tpins = (tgov + tcser + tparl + tpart)/4,
                      tptparl = (tpart + tparl)/2,
                      ybirth = ifelse(Q261 < 1, NA, Q261),
                      gencht = ifelse(ybirth < 1962, "indus",
                                      ifelse(ybirth %in% c(1962:1971), "386",
                                             ifelse(ybirth > 1971, "postd", NA)))) %>%
  select(strldr, strl2, jinsc, outs, aset, inc, age, gen, edu, ideo,
         apolr, region, hap, satfin, tpins, tptparl, gencht) %>% 
  mutate(strldrf = factor(strldr),
         region = factor(region, levels = c("Seoul", "Busan", "Daegu", "Incheon",
                                            "Gwangju", "Daejeon", "Ulsan",
                                            "Gyeonggi", "Gangwon", "Chungcheong",
                                            "Gyeongsang", "Jeolla")),
         regn = ifelse(region == "Seoul", 1,
                       ifelse(region == "Busan", 2,
                              ifelse(region == "Daegu", 3,
                                     ifelse(region == "Incheon", 4,
                                            ifelse(region == "Gwangju", 5,
                                                   ifelse(region == "Daejeon", 6,
                                                          ifelse(region == "Ulsan", 7,
                                                                 ifelse(region == "Gyeonggi", 8,
                                                                        ifelse(region == "Gangwon", 9,
                                                                               ifelse(region == "Chungcheong", 10,
                                                                                      ifelse(region == "Gyeongsang", 11, 12))))))))))))

library(MASS)
hp_1 <- polr(strldrf ~ jinsc + inc + age + gen + edu + ideo + apolr, data = wvs, Hess = TRUE)
hp_2 <- polr(strldrf ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)
hp_3 <- polr(strldrf ~ jinsc + I(inc^2) + inc + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)
hp_4 <- polr(strldrf ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)
hp_5 <- polr(strldrf ~ inc*jinsc + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)
hp_6 <- polr(strldrf ~ aset*outs + inc + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)
hp_7 <- polr(strldrf ~ inc*outs + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)

hl_1 <- lm(strldr ~ jinsc + inc + age + gen + edu + ideo + apolr, data = wvs)
hl_2 <- lm(strldr ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs)
hl_3 <- lm(strldr ~ jinsc + I(inc^2) + inc + age + gen + edu + ideo + apolr + region, data = wvs)
hl_4 <- lm(strldr ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs)
hl_5 <- lm(strldr ~ inc*jinsc + age + gen + edu + ideo + apolr + region, data = wvs)
hl_6 <- lm(strldr ~ aset*outs + inc + age + gen + edu + ideo + apolr + region, data = wvs)
hl_7 <- lm(strldr ~ inc*outs + age + gen + edu + ideo + apolr + region, data = wvs)

hg_1 <- glm(strl2 ~ jinsc + inc + age + gen + edu + ideo + apolr, data = wvs, family = binomial)
hg_2 <- glm(strl2 ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs, family = binomial)
hg_3 <- glm(strl2 ~ jinsc + I(inc^2) + inc + age + gen + edu + ideo + apolr + region, data = wvs, family = binomial)
hg_4 <- glm(strl2 ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs, family = binomial)
hg_5 <- glm(strl2 ~ inc*jinsc + age + gen + edu + ideo + apolr + region, data = wvs, family = binomial)
hg_6 <- glm(strl2 ~ aset*outs + inc + age + gen + edu + ideo + apolr + region, data = wvs, family = binomial)
hg_7 <- glm(strl2 ~ inc*outs + age + gen + edu + ideo + apolr + region, data = wvs, family = binomial)


# Figure 2

load("~/Desktop/Replication/WVS_TimeSeries_4_0.rdata")

wvt <- data1 %>% filter(COUNTRY_ALPHA == "KOR" & S002VS %in% c(3:7)) %>%
  select(E114, S020, S002VS, H006_01) %>% 
  mutate(strl = 5 - E114,
         strld = ifelse(strl == 1, "Very Bad",
                        ifelse(strl == 2, "Fairly Bad",
                               ifelse(strl == 3, "Fairly Good",
                                      ifelse(strl == 4, "Very Good)", NA)))),
         jins = 5 - H006_01) %>%
  rename(wave = S020)
 
f2a <- wvt %>% na.omit() %>% count(wave, strld) %>% group_by(wave) %>%
  mutate(r = n/sum(n)*100,
         strld = factor(strld, levels = c("Very Good)", "Fairly Good",
                                          "Fairly Bad", "Very Bad")),
         wave = factor(wave, levels = c("2018", "2010", "2005", "2001", "1996")))

f2at <- f2a %>% mutate(strl = ifelse(strld %in% c("Fairly Good", "Very Good)"),
                                     "Good", "Bad")) %>% 
  group_by(wave) %>% summarise(Bad = sum(n[strl == "Bad"]),
                               Good = sum(n[strl == "Good"])) %>% 
  mutate(r2 = round(Good/(Bad+Good)*100, 1))

library(ggplot2)
Fig2A <- ggplot(f2a, aes(x = wave, y = r, fill = strld)) + geom_col(width = 0.977) +
  coord_flip() +
  scale_fill_manual(values = c("#C8C9C7", "#A7A8AA", "#75787B", "#53565A")) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100),
                     labels = c("0%", "20%", "40%", "60%", "80%", "100%"),
                     expand = c(0,0)) +
  scale_x_discrete(expand = c(0,0)) +
  labs(fill = "(Strong Leader:  ") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  theme(axis.title.x = element_text(size = 5.5, face = 'bold', color = "white"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 11, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 11, face = 'bold', color = "black"),
        axis.ticks.y = element_blank()) +
  theme(legend.position = "top", legend.justification = "right",
        legend.title = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.text = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.margin = margin(t = 0, r = 0, b = -0.1, l = 0, unit = "cm"),
        plot.margin = margin(1, 15, 0, 5, "pt")) +
  guides(fill = guide_legend(reverse = TRUE))

for(i in 1:5){
  Fig2A <- Fig2A + annotate("text", x = f2at[[i, "wave"]],
                            y = 100 - f2at[[i, "r2"]],
                            label = paste0(f2at[[i, "r2"]], "%"),
                        hjust = -0.22, size = 3.5, color = "black")
  }

newd <- data.frame(jinsc = c(1:4), inc = c(rep(mean(wvs$inc, na.rm = T), 4)), 
                   age = c(rep(mean(wvs$age, na.rm = T), 4)),
                   gen = c(rep(mean(wvs$gen, na.rm = T), 4)),
                   edu = c(rep(mean(wvs$edu, na.rm = T), 4)),
                   ideo = c(rep(mean(wvs$ideo, na.rm = T), 4)),
                   apolr = c(rep(mean(wvs$apolr, na.rm = T), 4)))

predp <- predict(hg_1, newdata = newd, type = "response", se.fit = TRUE)

newd$predp <- predp$fit
newd$se <- predp$se.fit
newd$jinsc <- as.factor(newd$jinsc)

f2b <- wvt %>% filter(wave == 2018)

f2b <- f2b %>% select(jins, strld) %>% na.omit() %>% 
  mutate(strld = factor(strld, levels = c("Very Good)", "Fairly Good",
                                          "Fairly Bad", "Very Bad")))

f2b <- f2b %>% count(jins, strld) %>% group_by(jins) %>% mutate(r = n/sum(n)*100) %>% 
  mutate(jinsc = factor(jins))

f2b <- merge(f2b, newd, by = "jinsc") %>% mutate(jinsc = factor(jinsc))

Fig2B <- ggplot(f2b) +
  geom_col(aes(x = jinsc, y = r, fill = strld), width = 0.95, alpha = .8) +
  geom_point(aes(x = jinsc, y = predp*100), size = 1.8, shape = 15) +
  geom_errorbar(aes(x = jinsc, ymin = (predp - 1.96*se)*100,
                    ymax = (predp + 1.96*se)*100), width = 0.14) +
  scale_fill_manual(values = c("#C8C9C7", "#A7A8AA", "#75787B", "#53565A")) +
  scale_x_discrete(breaks = c(1:4), labels = c("1", "2", "3", "4"),
                   expand = c(0,0)) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100),
                     labels = c("0%", "20%", "40%", "60%", "80%", "100%"),
                     expand = c(0,0)) +
  labs(title = "2018", x = "Job Insecurity") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  theme(axis.title.x = element_text(size = 11, face = 'bold', color = "black"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 11, face = 'bold', color = "black",
                                   margin = margin(t = 0)),
        axis.text.y = element_text(size = 11, face = 'bold', color = "black"),
        axis.ticks.x = element_blank()) +
  theme(legend.position = "none",
        plot.margin = margin(1, 12, 3, 5, "pt"),
        plot.title = element_text(size = 11, face = 'bold', color = "black",
                                  hjust = 0.5, margin = margin(b = 11)))

library(gridExtra)
library(grid)

Fig2 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.025, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig2A,
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.1, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig2B,
    heights = c(0.1, 0.64)
  ),
  ncol = 2,
  widths = c(4, 1)
)

plot.new()
grid.draw(Fig2)


# Figure 3

coef <- coef(hp_4)
se <- sqrt(diag(vcov(hp_4)))
se <- head(se, length(se) - 3)
lb <- coef - 1.96*se
ub <- coef + 1.96*se

f3a <- data.frame(pred = names(coef), coef = coef, se = se, lb = lb, ub = ub) %>%
  filter(pred %in% c("aset", "jinsc", "inc", "age", "gen", "edu", "ideo", "apolr", "aset:jinsc")) %>%
  mutate(pred = ifelse(pred == "aset", "Assets",
                       ifelse(pred == "jinsc", "Job insecurity",
                              ifelse(pred == "inc", "Income",
                                     ifelse(pred == "age", "Age",
                                            ifelse(pred == "gen", "Gender",
                                                   ifelse(pred == "edu", "Education",
                                                          ifelse(pred == "ideo", "Ideological \n orientation",
                                                                 ifelse(pred == "apolr", "Affective \n polarization",
                                                                        "Assets × \n Job insecurity")))))))),
         pred = factor(pred, levels = rev(c("Job insecurity", "Assets",
                                            "Assets × \n Job insecurity", "Age",
                                            "Gender", "Education", "Income",
                                            "Ideological \n orientation",
                                            "Affective \n polarization"))),
         col = ifelse(lb*ub > 0, "#DA291C", "black")) %>% 
  arrange(pred)

Fig3A <- ggplot(f3a, aes(x = pred, y = coef, color = pred)) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", lwd = .2) +
  geom_point(size = 1.8, shape = 15) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.14) +
  coord_flip() +
  labs(x = NULL, y = "Coefficient Estimates") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.title.y = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.text.x = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 9.5, face = 'bold', color = "black"),
        legend.position = "none") +
  scale_color_manual(values = f3a$col)

coef <- coef(hp_5)
se <- sqrt(diag(vcov(hp_5)))
se <- head(se, length(se) - 3)
lb <- coef - 1.96*se
ub <- coef + 1.96*se

f3b <- data.frame(pred = names(coef), coef = coef, se = se, lb = lb, ub = ub) %>%
  filter(pred %in% c("inc", "jinsc", "age", "gen", "edu", "ideo", "apolr", "inc:jinsc")) %>%
  mutate(pred = ifelse(pred == "inc", "Income",
                       ifelse(pred == "jinsc", "Job insecurity",
                              ifelse(pred == "age", "Age",
                                     ifelse(pred == "gen", "Gender",
                                            ifelse(pred == "edu", "Education",
                                                   ifelse(pred == "ideo", "Ideological \n orientation",
                                                          ifelse(pred == "apolr", "Affective \n polarization",
                                                                 "Income × \n Job insecurity"))))))),
         pred = factor(pred, levels = rev(c("Job insecurity", "Income",
                                            "Income × \n Job insecurity", "Age",
                                            "Gender", "Education",
                                            "Ideological \n orientation",
                                            "Affective \n polarization"))),
         col = ifelse(lb*ub > 0, "#DA291C", "black")) %>% 
  arrange(pred)

Fig3B <- ggplot(f3b, aes(x = pred, y = coef, color = pred)) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", lwd = .2) +
  geom_point(size = 1.8, shape = 15) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.14) +
  coord_flip() +
  labs(x = NULL, y = "Coefficient Estimates") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.title.y = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.text.x = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 9.5, face = 'bold', color = "black"),
        legend.position = "none") +
  scale_color_manual(values = f3b$col)

Fig3 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.05, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig3A,
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.05, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig3B,
    heights = c(0.1, 0.75)
  ),
  ncol = 2,
  widths = c(1, 1)
)

plot.new()
grid.draw(Fig3)


# Figure 4

f4 <- as.matrix(wvs)
f4 <- as.data.frame(f4)
nc <- c("strl2", "jinsc", "aset", "inc", "age", "gen", "edu", "ideo", "apolr")
f4[, nc] <- apply(f4[nc], 2, as.numeric)

library(interflex)
out <- interflex(estimator = 'linear', method = 'logit',
                 Y = "strl2", D = "jinsc", X = "aset",
                 Z = c("age", "gen", "inc", "edu", "ideo", "apolr"),
                 data = f4, vcov.type = "robust", base = 1, na.rm = TRUE)

Fig4A <- plot(out, xlab = "Moderator: Assets",
              ylab = "Marginal Effect of Job Insecurity \n on Support for Undemocratic Leader \n",
              Xdistr = "density",
              cex.axis = 0.8, cex.lab = 0.8, theme.bw = TRUE, show.grid = FALSE, 
              order = c(2,3,4), subtitles = c("(Base: 1)", "2", "3", "4"),
              pool = TRUE, color = c("lightgray", "#A7A8AA", "black", "#DA291C"),
              cex.sub = 1, legend.title = "Job insecurity", show.subtitles = FALSE)

out <- interflex(estimator = 'linear', method = 'logit',
                 Y = "strl2", D = "jinsc", X = "inc",
                 Z = c("age", "gen", "edu", "ideo", "apolr"),
                 data = f4, vcov.type = "robust", base = 1, na.rm = TRUE)

Fig4B <- plot(out, xlab = "Moderator: Income", ylab = "", Xdistr = "density",
              cex.axis = 0.8, cex.lab = 0.8, theme.bw = TRUE, show.grid = FALSE,
              order = c(2,3,4), subtitles = c("(Base: 1)", "2", "3", "4"),
              pool = TRUE, color = c("lightgray", "#A7A8AA", "black", "#DA291C"),
              cex.sub = 1, legend.title = "Job insecurity")

Fig4 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.32, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig4A,
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.17, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig4B,
    heights = c(0.1, 0.75)
  ),
  ncol = 2, 
  widths = c(1, 1.23)
)

Fig4


# Table A5

ta5 <- wvs %>% select(-region, -strldrf, -strl2, -gencht, -regn)

dstat <- function(v) {
  n <- sum(!is.na(v))
  mean <- mean(v, na.rm = TRUE)
  sd <- sd(v, na.rm = TRUE)
  min <- min(v, na.rm = TRUE)
  med <- median(v, na.rm = TRUE)
  max <- max(v, na.rm = TRUE)
  
  des <- data.frame(
    n = n,
    mean = mean,
    sd = sd,
    min = as.character(min),
    median = as.character(med),
    max = as.character(max)
  )
  
  return(des)
}

dslist <- lapply(ta5, dstat)
dstat <- do.call(rbind, dslist)
dstat <- rownames_to_column(dstat, var = 'variable')

library(xtable)
TabA5 <- xtable(dstat, digits = 3)

table(wvs$region)

TabA5


# Table A7

round(cor(wvs$jinsc, wvs$aset, use = "complete.obs"), 3)
round(cor(wvs$jinsc, wvs$inc, use = "complete.obs"), 3)

ta7 <- wvs %>% group_by(jinsc, aset) %>% summarise(n = n(), r = round(n/1245*100, 1)) %>% print(n = 33)

ta7 %>% group_by(jinsc) %>% summarise(n = sum(n), r = round(n/1245*100, 1))

ta7 %>% group_by(aset) %>% summarise(n = sum(n), r = round(n/1245*100, 1))


# Table A8

ta8 <- wvs %>% group_by(jinsc, inc) %>% summarise(n = n(), r = round(n/1245*100, 1)) %>% print(n = 33)

ta8 %>% group_by(jinsc) %>% summarise(n = sum(n), r = round(n/1245*100, 1))

ta8 %>% group_by(inc) %>% summarise(n = sum(n), r = round(n/1245*100, 1))


# Table A10

library(stargazer)
stargazer(hp_1, hp_2, hp_3, hp_4, hp_5, hp_6, hp_7, type = "latex",
          title = "Job insecurity, assets, income, and undemocratic attitudes (WVS)",
          star.char = c("+", "*", "**"))


# Table A11

stargazer(hl_1, hl_2, hl_3, hl_4, hl_5, hl_6, hl_7, type = "latex",
          title = "Job insecurity, assets, income, and undemocratic attitudes (WVS, OLS)",
          star.char = c("+", "*", "**"))


# Table A12

stargazer(hg_1, hg_2, hg_3, hg_4, hg_5, hg_6, hg_7, type = "latex",
          title = "Job insecurity, assets, income, and undemocratic attitudes (WVS, GLM)",
          star.char = c("+", "*", "**"))


# Table A14

ta14_1 <- polr(as.factor(hap) ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs, Hess = TRUE)

ta14_2 <- lm(satfin ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs)

ta14_3 <- lm(tpins ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs)

ta14_4 <- lm(tptparl ~ jinsc + inc + age + gen + edu + ideo + apolr + region, data = wvs)

stargazer(ta14_1, ta14_2, ta14_3, ta14_4, type = "latex",
          title = "Mechanism test: consequences of job insecurity (WVS)",
          star.char = c("+", "*", "**"))


# Table A15

ta15_1 <- polr(strldrf ~ jinsc + inc + age + gen + edu + ideo + apolr + region + gencht,
               data = wvs, Hess = TRUE)

ta15_2 <- polr(strldrf ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + region + gencht,
               data = wvs, Hess = TRUE)

ta15_3 <- polr(strldrf ~ inc*jinsc + age + gen + edu + ideo + apolr + region + gencht,
               data = wvs, Hess = TRUE)

stargazer(ta15_1, ta15_2, ta15_3, type = "latex",
          title = "Additional control: generational cohort dummies",
          star.char = c("+", "*", "**"))


# Figure A1

fa123 <- data1 %>% filter(COUNTRY_ALPHA == "KOR") %>% filter(S002VS %in% c(3:7)) %>%
  select(E114, S020, E179_WVS7LOC, X001, X003) %>% 
  mutate(strl = 5 - E114,
         strld = ifelse(strl == 1, "Very Bad",
                        ifelse(strl == 2, "Fairly Bad",
                               ifelse(strl == 3, "Fairly Good",
                                      ifelse(strl == 4, "Very Good)", NA)))),
         prty = ifelse(E179_WVS7LOC == 410001, "Conservative",
                       ifelse(E179_WVS7LOC == 410002, "Liberal", NA)),
         gen = ifelse(X001 == 1, "Male", ifelse(X001 == 2, "Female", NA)),
         ageg = ifelse(X003 %in% c(20:39), "20s-30s",
                       ifelse(X003 %in% c(40:59), "40s-50s",
                              ifelse(X003 %in% c(60:79), "60s-70s", NA)))) %>% 
  rename(wave = S020)

fa1 <- fa123 %>% filter(wave %in% c(2005:2018))

fa1_1 <- fa1 %>% filter(prty == "Conservative")

fa1_1 <- fa1_1 %>% na.omit() %>% count(wave, strld) %>% group_by(wave) %>%
  mutate(r = n/sum(n)*100,
         strld = factor(strld, levels = c("Very Good)", "Fairly Good",
                                          "Fairly Bad", "Very Bad")),
         wave = factor(wave, levels = c("2018", "2010", "2005")))

fa1_1t <- fa1_1 %>% mutate(strl = ifelse(strld %in% c("Fairly Good", "Very Good)"),
                                     "Good", "Bad")) %>% 
  group_by(wave) %>% summarise(Bad = sum(n[strl == "Bad"]),
                               Good = sum(n[strl == "Good"])) %>% 
  mutate(r2 = round(Good/(Bad+Good)*100, 1))

FigA1a <- ggplot(fa1_1, aes(x = wave, y = r, fill = strld)) + geom_col(width = 0.977) +
  coord_flip() +
  scale_fill_manual(values = c("#C8C9C7", "#A7A8AA", "#75787B", "#53565A")) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100),
                     labels = c("0%", "20%", "40%", "60%", "80%", "100%"),
                     expand = c(0,0)) +
  scale_x_discrete(expand = c(0,0)) +
  labs(fill = "(Strong Leader:  ", title = "Conservative") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  theme(axis.title.x = element_text(size = 5.5, face = 'bold', color = "white"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 11, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 11, face = 'bold', color = "black"),
        axis.ticks.y = element_blank()) +
  theme(legend.position = "top", legend.justification = "right",
        legend.title = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.text = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.margin = margin(t = -0.8, r = 0, b = -0.1, l = 0, unit = "cm"),
        plot.margin = margin(5, 15, 0, 5, "pt"),
        plot.title = element_text(size = 14, face = 'bold', color = "black")) +
  guides(fill = guide_legend(reverse = TRUE))

for(i in 1:3){
  FigA1a <- FigA1a + annotate("text", x = fa1_1t[[i, "wave"]],
                            y = 100 - fa1_1t[[i, "r2"]],
                            label = paste0(fa1_1t[[i, "r2"]], "%"),
                            hjust = -0.22, size = 3.5, color = "black")
}

FigA1a

fa1_2 <- fa1 %>% filter(prty == "Liberal")

fa1_2 <- fa1_2 %>% na.omit() %>% count(wave, strld) %>% group_by(wave) %>%
  mutate(r = n/sum(n)*100,
         strld = factor(strld, levels = c("Very Good)", "Fairly Good",
                                          "Fairly Bad", "Very Bad")),
         wave = factor(wave, levels = c("2018", "2010", "2005")))

fa1_2t <- fa1_2 %>% mutate(strl = ifelse(strld %in% c("Fairly Good", "Very Good)"),
                                         "Good", "Bad")) %>% 
  group_by(wave) %>% summarise(Bad = sum(n[strl == "Bad"]),
                               Good = sum(n[strl == "Good"])) %>% 
  mutate(r2 = round(Good/(Bad+Good)*100, 1))

FigA1b <- ggplot(fa1_2, aes(x = wave, y = r, fill = strld)) + geom_col(width = 0.977) +
  coord_flip() +
  scale_fill_manual(values = c("#C8C9C7", "#A7A8AA", "#75787B", "#53565A")) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100),
                     labels = c("0%", "20%", "40%", "60%", "80%", "100%"),
                     expand = c(0,0)) +
  scale_x_discrete(expand = c(0,0)) +
  labs(fill = "(Strong Leader:  ", title = "Liberal") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  theme(axis.title.x = element_text(size = 5.5, face = 'bold', color = "white"),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 11, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 11, face = 'bold', color = "black"),
        axis.ticks.y = element_blank()) +
  theme(legend.position = "top", legend.justification = "right",
        legend.title = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.text = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.margin = margin(t = -0.8, r = 0, b = -0.1, l = 0, unit = "cm"),
        plot.margin = margin(5, 15, 0, 5, "pt"),
        plot.title = element_text(size = 14, face = 'bold', color = "black")) +
  guides(fill = guide_legend(reverse = TRUE))

for(i in 1:3){
  FigA1b <- FigA1b + annotate("text", x = fa1_2t[[i, "wave"]],
                              y = 100 - fa1_2t[[i, "r2"]],
                              label = paste0(fa1_2t[[i, "r2"]], "%"),
                              hjust = -0.22, size = 3.5, color = "black")
}

FigA1b


# Figure A2

fa23 <- fa123 %>% select(-prty)

disagg <- function(v, vt){
  f <- fa23 %>% filter({{v}} == vt)
  
  f <- f %>% na.omit() %>%  count(wave, strld) %>% group_by(wave) %>%
    mutate(r = n/sum(n)*100,
           strld = factor(strld, levels = c("Very Good)", "Fairly Good",
                                            "Fairly Bad", "Very Bad")),
           wave = factor(wave, levels = c("2018", "2010", "2005", "2001", "1996")))
  
  ft <- f %>% mutate(strl = ifelse(strld %in% c("Fairly Good", "Very Good)"),
                                           "Good", "Bad")) %>% 
    group_by(wave) %>% summarise(Bad = sum(n[strl == "Bad"]),
                                 Good = sum(n[strl == "Good"])) %>% 
    mutate(r2 = round(Good/(Bad+Good)*100, 1))
  
  fig <- ggplot(f, aes(x = wave, y = r, fill = strld)) + geom_col(width = 0.977) +
    coord_flip() +
    scale_fill_manual(values = c("#C8C9C7", "#A7A8AA", "#75787B", "#53565A")) +
    scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100),
                       labels = c("0%", "20%", "40%", "60%", "80%", "100%"),
                       expand = c(0,0)) +
    scale_x_discrete(expand = c(0,0)) +
    labs(fill = "(Strong Leader:  ", title = vt) +
    theme_bw() +
    theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
    theme(axis.title.x = element_text(size = 5.5, face = 'bold', color = "white"),
          axis.title.y = element_blank(),
          axis.text.x = element_text(size = 11, face = 'bold', color = "black"),
          axis.text.y = element_text(size = 11, face = 'bold', color = "black"),
          axis.ticks.y = element_blank()) +
    theme(legend.position = "top", legend.justification = "right",
          legend.title = element_text(size = 10.5, face = 'bold', color = "black"),
          legend.text = element_text(size = 10.5, face = 'bold', color = "black"),
          legend.margin = margin(t = -0.8, r = 0, b = -0.1, l = 0, unit = "cm"),
          plot.margin = margin(5, 15, 0, 5, "pt"),
          plot.title = element_text(size = 14, face = 'bold', color = "black")) +
    guides(fill = guide_legend(reverse = TRUE))
  
  for(i in 1:5){
    fig <- fig + annotate("text", x = ft[[i, "wave"]],
                          y = 100 - ft[[i, "r2"]],
                          label = paste0(ft[[i, "r2"]], "%"),
                          hjust = -0.22, size = 3.5, color = "black")
  }
  
  return(fig)
}

FigA2a <- disagg(gen, "Male")
FigA2b <- disagg(gen, "Female")

FigA2a
FigA2b


# Figure A3

FigA3a <- disagg(ageg, "20s-30s")
FigA3b <- disagg(ageg, "40s-50s")
FigA3c <- disagg(ageg, "60s-70s")

FigA3a
FigA3b
FigA3c


# Figure A4

jknifL <- function(dat, v, vt, rmod, mtyp, var, varn, ordr){
  dat <- dat %>% filter({{v}} != vt)
  
  if (mtyp == "lm") {
    jk <- lm(rmod, data = dat)
  } else if (mtyp == "glm") {
    jk <- glm(rmod, data = dat, family = binomial)
  } else if (mtyp == "polr") {
    require(MASS)
    jk <- polr(rmod, data = dat, Hess = TRUE)
  } else {
    stop("Support 'lm', 'glm', 'polr'")
  }
  
  coef <- coef(jk)
  se <- sqrt(diag(vcov(jk)))
  
  if (mtyp %in% c("lm", "glm")) {
    se <- head(se, length(se))
  } else if (mtyp == "polr") {
    se <- head(se, length(se) - 3)
  }
  
  lb <- coef - 1.96*se
  ub <- coef + 1.96*se
  
  jd <- data.frame(pred = names(coef), coef = coef, se = se, lb = lb, ub = ub) %>%
    filter(pred %in% var)
  
  nam <- data.frame(var, varn) %>% rename(pred = var)
  
  jd <- merge(jd, nam, by = "pred") %>% select(-pred) %>% 
    mutate(varn = factor(varn, levels = rev(ordr)),
           col = ifelse(lb*ub > 0, "#DA291C", "black")) %>% 
    arrange(varn)
  
  plot <- ggplot(jd, aes(x = varn, y = coef, color = varn)) +
    geom_hline(yintercept = 0, linetype = "solid", color = "black", lwd = .2) +
    geom_point(size = 1.8, shape = 15) +
    geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.14) +
    coord_flip() +
    labs(x = NULL, y = "Coefficient Estimates",
         title = paste0("(Without ", vt, ")")) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 9, face = 'bold', color = "black"),
          axis.title.y = element_text(size = 9, face = 'bold', color = "black"),
          axis.text.x = element_text(size = 9, face = 'bold', color = "black"),
          axis.text.y = element_text(size = 9, face = 'bold', color = "black"),
          legend.position = "none",
          plot.title = element_text(hjust = 1, size = 9, face = 'bold', color = "black")) +
    scale_color_manual(values = jd$col)
  
  return(plot)
}

jknifR <- function(dat, v, vt, rmod, mtyp, var, varn, ordr){
  dat <- dat %>% filter({{v}} != vt)
  
  if (mtyp == "lm") {
    jk <- lm(rmod, data = dat)
  } else if (mtyp == "glm") {
    jk <- glm(rmod, data = dat, family = binomial)
  } else if (mtyp == "polr") {
    require(MASS)
    jk <- polr(rmod, data = dat, Hess = TRUE)
  } else {
    stop("Support 'lm', 'glm', 'polr'")
  }
  
  coef <- coef(jk)
  se <- sqrt(diag(vcov(jk)))
  
  if (mtyp %in% c("lm", "glm")) {
    se <- head(se, length(se))
  } else if (mtyp == "polr") {
    se <- head(se, length(se) - 3)
  }
  
  lb <- coef - 1.96*se
  ub <- coef + 1.96*se
  
  jd <- data.frame(pred = names(coef), coef = coef, se = se, lb = lb, ub = ub) %>%
    filter(pred %in% var)
  
  nam <- data.frame(var, varn) %>% rename(pred = var)
  
  jd <- merge(jd, nam, by = "pred") %>% select(-pred) %>% 
    mutate(varn = factor(varn, levels = rev(ordr)),
           col = ifelse(lb*ub > 0, "#DA291C", "black")) %>% 
    arrange(varn)
  
  plot <- ggplot(jd, aes(x = varn, y = coef, color = varn)) +
    geom_hline(yintercept = 0, linetype = "solid", color = "black", lwd = .2) +
    geom_point(size = 1.8, shape = 15) +
    geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.14) +
    coord_flip() +
    labs(x = NULL, y = "Coefficient Estimates",
         title = paste0("(Without ", vt, ")")) +
    theme_bw() +
    theme(axis.title.x = element_text(size = 9, face = 'bold', color = "black"),
          axis.title.y = element_text(size = 9, face = 'bold', color = "black"),
          axis.text.x = element_text(size = 9, face = 'bold', color = "black"),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          legend.position = "none",
          plot.title = element_text(hjust = 1, size = 9, face = 'bold', color = "black")) +
    scale_color_manual(values = jd$col)
  
  return(plot)
}

regions <- names(table(wvs$region))
plist <- list()

for(i in c(seq(1, 10, by = 3))){
  plt <- jknifL(wvs, region, regions[i],
                strldrf ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + as.factor(regn),
                "polr",
                c("aset", "jinsc", "inc", "age", "gen", "edu", "ideo", "apolr", "aset:jinsc"),
                c("Assets", "Job insecurity", "Income", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Assets × Job insecurity"),
                c("Job insecurity", "Assets", "Assets × Job insecurity", "Age",
                  "Gender", "Education", "Income", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

for(i in c(2,3,5,6,8,9,11,12)){
  plt <- jknifR(wvs, region, regions[i],
                strldrf ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + as.factor(regn),
                "polr",
                c("aset", "jinsc", "inc", "age", "gen", "edu", "ideo", "apolr", "aset:jinsc"),
                c("Assets", "Job insecurity", "Income", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Assets × Job insecurity"),
                c("Job insecurity", "Assets", "Assets × Job insecurity", "Age",
                  "Gender", "Education", "Income", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

FigA4 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[1]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[2]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("C", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[3]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("D", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[4]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("E", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[5]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("F", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[6]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("G", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[7]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("H", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[8]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("I", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[9]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("J", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[10]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("K", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[11]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("L", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[12]],
    heights = c(0.1, 0.75)
  ),
  ncol = 3,
  widths = c(1.642, 1, 1)
)

plot.new()
grid.draw(FigA4)


# Figure A6

wvs$ageg <- ifelse(wvs$age == 1, "~ 24",
                   ifelse(wvs$age == 2, "25 ~ 34",
                          ifelse(wvs$age == 3, "35 ~ 44",
                                 ifelse(wvs$age == 4, "45 ~ 54",
                                        ifelse(wvs$age == 5, "55 ~ 64", "65 ~ ")))))

agegrp <- names(table(wvs$ageg))
plist <- list()

for(i in c(1,4)){
  plt <- jknifL(wvs, ageg, agegrp[i],
                strldrf ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + as.factor(regn),
                "polr",
                c("aset", "jinsc", "inc", "age", "gen", "edu", "ideo", "apolr", "aset:jinsc"),
                c("Assets", "Job insecurity", "Income", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Assets × Job insecurity"),
                c("Job insecurity", "Assets", "Assets × Job insecurity", "Age",
                  "Gender", "Education", "Income", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

for(i in c(2,3,5,6)){
  plt <- jknifR(wvs, ageg, agegrp[i],
                strldrf ~ aset*jinsc + inc + age + gen + edu + ideo + apolr + as.factor(regn),
                "polr",
                c("aset", "jinsc", "inc", "age", "gen", "edu", "ideo", "apolr", "aset:jinsc"),
                c("Assets", "Job insecurity", "Income", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Assets × Job insecurity"),
                c("Job insecurity", "Assets", "Assets × Job insecurity", "Age",
                  "Gender", "Education", "Income", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

FigA6 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[1]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[2]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("C", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[3]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("D", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[4]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("E", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[5]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("F", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[6]],
    heights = c(0.1, 0.75)
  ),
  ncol = 3,
  widths = c(1.642, 1, 1)
)

plot.new()
grid.draw(FigA6)


#####
# KGSS
# Figure 5
# Table A6, Table A9, Table A13, Table A15
# Figure A5, Figure A7

kgss <- read_sav("~/Desktop/Replication/2003-2021_KGSS_public_v3_20230210.sav")
select <- dplyr::select

kg <- kgss %>% filter(YEAR == 2021) %>%
  mutate(demdic = ifelse(DEMODICT == 1, 0, ifelse(DEMODICT %in% c(2:3), 1, NA)),
         jobinsc = ifelse(WGSTAB == 1, 0, ifelse(WGSTAB == 2, 1, NA)),
         WGSTAT = ifelse(WGSTAT < 1, NA, WGSTAT), WGPTFL = ifelse(WGPTFL < 1, NA, WGPTFL),
         lous = ifelse(WGSTAT %in% c(2,3), 1, 0), lous2 = ifelse(WGPTFL == 1, 1, 0),
         lous = ifelse(lous2 == 1, 1, lous),
         income = ifelse(INCOM0 < 0, NA, INCOM0),
         hompop = ifelse(HOMPOP < 1, NA, HOMPOP),
         adjinc = income / sqrt(hompop),
         inc = ifelse(adjinc < 753/12, 1,
                      ifelse(adjinc < 1541/12, 2,
                             ifelse(adjinc < 2123/12, 3,
                                    ifelse(adjinc < 2678/12, 4,
                                           ifelse(adjinc < 3204/12, 5,
                                                  ifelse(adjinc < 3816/12, 6,
                                                         ifelse(adjinc < 4537/12, 7,
                                                                ifelse(adjinc < 5555/12, 8,
                                                                       ifelse(adjinc < 7327/12, 9, 10))))))))),
         age = ifelse(AGE < 1, NA,
                      ifelse(AGE < 30, 1,
                             ifelse(AGE < 40, 2,
                                    ifelse(AGE < 50, 3,
                                           ifelse(AGE < 60, 4,
                                                  ifelse(AGE < 70, 5, 6)))))),
         gen = ifelse(SEX == 2, 0, SEX),
         edu = ifelse(EDUC < 0 | EDUC > 7, NA, EDUC),
         ideo = ifelse(PARTYLR < 1, NA, PARTYLR),
         fmj = ifelse(LIKEPRTY1 < 0, NA, LIKEPRTY1),
         fkh = ifelse(LIKEPRTY2 < 0, NA, LIKEPRTY2),
         apolr = abs(fmj - fkh),
         regn = REGION,
         region = ifelse(REGION == 1, "Seoul",
                         ifelse(REGION == 2, "Gyeonggi",
                                ifelse(REGION == 3, "Gangwon",
                                       ifelse(REGION == 4, "Chungcheong",
                                              ifelse(REGION == 5, "Gyeongsang",
                                                     ifelse(REGION == 6, "Jeolla",
                                                            "Jeju")))))),
         gencht = ifelse(AGE < 51, "postd",
                         ifelse(AGE %in% c(51:60), "386",
                                ifelse(AGE > 60, "indus", NA)))) %>% 
  select(demdic, jobinsc, lous, inc, age, gen, edu, ideo, apolr,
         region, regn, gencht) %>% 
  mutate(region = factor(region, levels = c("Seoul", "Gyeonggi", "Gangwon",
                                            "Chungcheong", "Gyeongsang", "Jeolla", "Jeju")))

hg1 <- glm(demdic ~ jobinsc + inc + age + gen + edu + ideo + apolr, data = kg, family = binomial)
hg2 <- glm(demdic ~ jobinsc + inc + age + gen + edu + ideo + apolr + region, data = kg, family = binomial)
hg3 <- glm(demdic ~ jobinsc + I(inc^2) + inc + age + gen + edu + ideo + apolr + region, data = kg, family = binomial)
hg4 <- glm(demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + region, data = kg, family = binomial)
hg5 <- glm(demdic ~ inc*lous + age + gen + edu + ideo + apolr + region, data = kg, family = binomial)

hl1 <- lm(demdic ~ jobinsc + inc + age + gen + edu + ideo + apolr, data = kg)
hl2 <- lm(demdic ~ jobinsc + inc + age + gen + edu + ideo + apolr + region, data = kg)
hl3 <- lm(demdic ~ jobinsc + I(inc^2) + inc + age + gen + edu + ideo + apolr + region, data = kg)
hl4 <- lm(demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + region, data = kg)
hl5 <- lm(demdic ~ inc*lous + age + gen + edu + ideo + apolr + region, data = kg)


# Figure 5

coef <- coef(hg4)
se <- sqrt(diag(vcov(hg4)))
lb <- coef - 1.96*se
ub <- coef + 1.96*se

f5a <- data.frame(pred = names(coef), coef = coef, se = se, lb = lb, ub = ub) %>% 
  filter(pred %in% c("inc", "jobinsc", "age", "gen", "edu", "ideo", "apolr", "inc:jobinsc")) %>% 
  mutate(pred = ifelse(pred == "inc", "Income",
                       ifelse(pred == "jobinsc", "Job insecurity",
                              ifelse(pred == "age", "Age",
                                     ifelse(pred == "gen", "Gender",
                                            ifelse(pred == "edu", "Education",
                                                   ifelse(pred == "ideo", "Ideological \n orientation",
                                                          ifelse(pred == "apolr", "Affective \n polarization",
                                                                 "Income × \n Job insecurity"))))))),
         pred = factor(pred, levels = rev(c("Job insecurity", "Income",
                                            "Income × \n Job insecurity", "Age",
                                            "Gender", "Education",
                                            "Ideological \n orientation",
                                            "Affective \n polarization"))),
         col = ifelse(lb*ub > 0, "#DA291C", "black")) %>% 
  arrange(pred)
  
Fig5A <- ggplot(f5a, aes(x = pred, y = coef, color = pred)) +
  geom_hline(yintercept = 0, linetype = "solid", color = "black", lwd = .2) +
  geom_point(size = 1.8, shape = 15) +
  geom_errorbar(aes(ymin = lb, ymax = ub), width = 0.14) +
  coord_flip() +
  labs(x = NULL, y = "Coefficient Estimates") +
  theme_bw() +
  theme(axis.title.x = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.title.y = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.text.x = element_text(size = 9.5, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 9.5, face = 'bold', color = "black"),
        legend.position = "none") +
  scale_color_manual(values = f5a$col)

f5 <- as.matrix(kg)
f5 <- as.data.frame(f5)
nc <- c("demdic", "jobinsc", "inc", "age", "gen", "edu", "ideo", "apolr")
f5[, nc] <- apply(f5[nc], 2, as.numeric)

out <- interflex(estimator = 'linear', method = 'logit',
                 Y = "demdic", D = "jobinsc", X = "inc",
                 Z = c("age", "gen", "edu", "ideo", "apolr"),
                 data = f5, vcov.type = "robust", base = 0, na.rm = TRUE)

Fig5B <- plot(out, xlab = "Moderator: Income",
              ylab = "Marginal Effect of Job Insecurity \n on Support for Undemocratic Regime",
              Xdistr = "density", cex.axis = 0.7, cex.lab = 0.7,
              theme.bw = TRUE, show.grid = FALSE, diff.values = c(1, 5, 10),
              order = c(1), subtitles = c("(Base: secure)", "Insecure"),
              pool = TRUE, color = c("#C8C9C7", "black"), cex.sub = 1,
              legend.title = "Treatment: Job")

empty_plot <- grid.draw(rectGrob(gp = gpar(fill = NA, col = NA)))

Fig5 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.06, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig5A,
    heights = c(0.1, 0.75)
  ),
  empty_plot,
  arrangeGrob(
    textGrob("B", x = 0.035, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    Fig5B,
    heights = c(0.1, 0.75)
  ),
  ncol = 3, 
  widths = c(1, 0.05, 1.25)
)

plot.new()
grid.draw(Fig5)


# Table A6

ta6 <- kg %>% select(-region, -regn, -gencht)

dstat <- function(v) {
  n <- sum(!is.na(v))
  mean <- mean(v, na.rm = TRUE)
  sd <- sd(v, na.rm = TRUE)
  min <- min(v, na.rm = TRUE)
  med <- median(v, na.rm = TRUE)
  max <- max(v, na.rm = TRUE)
  
  des <- data.frame(
    n = n,
    mean = mean,
    sd = sd,
    min = as.character(min),
    median = as.character(med),
    max = as.character(max)
  )
  
  return(des)
}

dslist <- lapply(ta6, dstat)
dstat <- do.call(rbind, dslist)
dstat <- rownames_to_column(dstat, var = 'variable')

library(xtable)
TabA6 <- xtable(dstat, digits = 3)

table(kg$region)

TabA6


# Table A9

ta9 <- kg %>% group_by(jobinsc, inc) %>% summarise(n = n(), r = round(n/1205 * 100, 1)) %>% print(n = 33)

ta9 %>% group_by(jobinsc) %>% summarise(num = sum(n), rat = round(num/1205*100, 1))

ta9 %>% group_by(inc) %>% summarise(num = sum(n), rat = round(num/1205*100, 1))


# Table A13

stargazer(hg1, hg2, hg3, hg4, hg5, hl1, hl2, hl3, hl4, hl5, type = "latex",
          title = "Job insecurity, income, and undemocratic attitudes (KGSS)",
          star.char = c("+", "*", "**"))


# Table A15

ta15_4 <- glm(demdic ~ jobinsc + inc + age + gen + edu + ideo + apolr + region + gencht,
              data = kg, family = binomial)

ta15_5 <- glm(demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + region + gencht,
              data = kg, family = binomial)

stargazer(ta15_4, ta15_5, type = "latex",
          title = "Additional control: generational cohort dummies",
          star.char = c("+", "*", "**"))


# Figure A5

regions <- names(table(kg$region))
plist <- list()

for(i in c(1, 4, 7)){
  plt <- jknifL(kg, region, regions[i],
                demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + as.factor(regn),
                "glm",
                c("inc", "jobinsc", "age", "gen", "edu", "ideo", "apolr", "inc:jobinsc"),
                c("Income", "Job insecurity", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Income × Job insecurity"),
                c("Job insecurity", "Income", "Income × Job insecurity", "Age",
                  "Gender", "Education", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

for(i in c(2, 3, 5, 6)){
  plt <- jknifR(kg, region, regions[i],
                demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + as.factor(regn),
                "glm",
                c("inc", "jobinsc", "age", "gen", "edu", "ideo", "apolr", "inc:jobinsc"),
                c("Income", "Job insecurity", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Income × Job insecurity"),
                c("Job insecurity", "Income", "Income × Job insecurity", "Age",
                  "Gender", "Education", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

FigA5 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[1]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[2]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("C", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[3]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("D", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[4]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("E", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[5]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("F", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[6]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("G", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[7]],
    heights = c(0.1, 0.75)
  ),
  ncol = 3,
  widths = c(1.642, 1, 1)
)

plot.new()
grid.draw(FigA5)


# Figure A7

kg$ageg <- ifelse(kg$age == 1, "~ 29",
                  ifelse(kg$age == 2, "30 ~ 39",
                         ifelse(kg$age == 3, "40 ~ 49",
                                ifelse(kg$age == 4, "50 ~ 59",
                                       ifelse(kg$age == 5, "60 ~ 69", "70 ~ ")))))

agegrp <- names(table(kg$ageg))
plist <- list()

for(i in c(1, 4)){
  plt <- jknifL(kg, ageg, agegrp[i],
                demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + as.factor(regn),
                "glm",
                c("inc", "jobinsc", "age", "gen", "edu", "ideo", "apolr", "inc:jobinsc"),
                c("Income", "Job insecurity", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Income × Job insecurity"),
                c("Job insecurity", "Income", "Income × Job insecurity", "Age",
                  "Gender", "Education", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

for(i in c(2, 3, 5, 6)){
  plt <- jknifR(kg, ageg, agegrp[i],
                demdic ~ inc*jobinsc + age + gen + edu + ideo + apolr + as.factor(regn),
                "glm",
                c("inc", "jobinsc", "age", "gen", "edu", "ideo", "apolr", "inc:jobinsc"),
                c("Income", "Job insecurity", "Age", "Gender", "Education",
                  "Ideological orientation", "Affective polarization", "Income × Job insecurity"),
                c("Job insecurity", "Income", "Income × Job insecurity", "Age",
                  "Gender", "Education", "Ideological orientation",
                  "Affective polarization"))
  plist[[i]] <- plt
}

FigA7 <- grid.arrange(
  arrangeGrob(
    textGrob("A", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[1]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("B", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[2]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("C", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[3]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("D", x = 0.435, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[4]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("E", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[5]],
    heights = c(0.1, 0.75)
  ),
  arrangeGrob(
    textGrob("F", x = 0.085, y = 0.8, just = "top", gp = gpar(fontsize = 16)),
    plist[[6]],
    heights = c(0.1, 0.75)
  ),
  ncol = 3,
  widths = c(1.642, 1, 1)
)

plot.new()
grid.draw(FigA7)


#####
# Figure A8

library(readr)
dic <- read_csv("~/Desktop/Replication/sentiment.csv")

m_comment <- read_csv("dp_cleaned.csv")

library(tidytext)
library(stringr)
word_comment <- m_comment %>%
  unnest_tokens(input = comment, output = word, token = "words", drop = F) %>% 
  filter(str_detect(word, "[가-힣]") & str_count(word) >= 2)

word_comment <- word_comment %>% left_join(dic, by = "word",
                                           relationship = "many-to-many") %>% 
  mutate(polarity = ifelse(is.na(polarity), 0, polarity))

score_comment <- word_comment %>% group_by(id, comment, year) %>%
  summarise(score = sum(polarity)) %>% ungroup()

score_comment <- score_comment %>% mutate(sentiment = ifelse(score >= 1, "pos",
                                                             ifelse(score <= -1, "neg",
                                                                    "neu")))

frequency_score <- score_comment %>% group_by(year) %>%
  count(sentiment) %>% mutate(ratio = n/sum(n)*100) %>% 
  mutate(party = "dp")

k_comment <- read_csv("ppp_cleaned.csv")

word_comment2 <- k_comment %>%
  unnest_tokens(input = comment, output = word, token = "words", drop = F) %>% 
  filter(str_detect(word, "[가-힣]") & str_count(word) >= 2)

word_comment2 <- word_comment2 %>% left_join(dic, by = "word") %>% 
  mutate(polarity = ifelse(is.na(polarity), 0, polarity))

score_comment2 <- word_comment2 %>% group_by(id, comment, year) %>%
  summarise(score = sum(polarity)) %>% ungroup()

score_comment2 <- score_comment2 %>% mutate(sentiment = ifelse(score >= 1, "pos",
                                                               ifelse(score <= -1,
                                                                      "neg", "neu")))

frequency_score2 <- score_comment2 %>% group_by(year) %>%
  count(sentiment) %>% mutate(ratio = n/sum(n)*100) %>% 
  mutate(party = "ppp")

score <- rbind(frequency_score, frequency_score2) %>% 
  mutate(yr = paste0("'", substr(year, 3,4)),
         party = ifelse(party == "dp", "Democratic Party", "People Power Party"),
         sentiment = ifelse(sentiment == "pos", "Positive",
                            ifelse(sentiment == "neu", "Neutral", "Negative)")),
         sentiment = factor(sentiment, levels = c("Positive", "Neutral", "Negative)")))

FigA8 <- ggplot(score, aes(x = yr, y = ratio, fill = sentiment)) +
  geom_col(width = 0.9475) +
  scale_fill_manual(values = c("#C8C9C7", "#A7A8AA", "#75787B")) +
  scale_x_discrete(expand = c(0,0)) +
  scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100),
                     labels = c("0%", "20%", "40%", "60%", "80%", "100%"),
                     expand = c(0,0)) +
  labs(fill = "(Comment Sentiments:  ") +
  theme_bw() +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank()) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 10, face = 'bold', color = "black"),
        axis.text.y = element_text(size = 10, face = 'bold', color = "black"),
        axis.ticks.x = element_blank()) +
  theme(panel.spacing.x = unit(0.45, "cm"),
        legend.position = "top", legend.justification = "right",
        legend.title = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.text = element_text(size = 10.5, face = 'bold', color = "black"),
        legend.margin = margin(t = 0, r = 0, b = -0.1, l = 0, unit = "cm")) +
  facet_grid(~ party) +
  theme(strip.text = element_text(size = 15, face = "bold", color = "black"),
        strip.background = element_rect(fill = "transparent", color = "transparent"))

FigA8

